#include "dyn_roots_hunting.h"

Hunting::Hunting(Functor& _func,
                 Functor& _derivative,
                 double _initialPoint,
                 double _endPoint,
                 double _scanningSpeed,
                 double _tolBis,
                 double _tolQuad,
                 double _maxValue,
                 int _maxIterations)
    : function(_func), derivative(_derivative)
{
    setInitialValue(_initialPoint);
    setEndPoint(_endPoint);
    setScanningSpeed(_scanningSpeed);
    setToleranceBis(_tolBis);
    setTolerance(_tolQuad);
    setMaxValue(_maxValue);
    setMaxIterations(_maxIterations);
}

int Hunting::solve()
{
    // We check to see if b>a. If not, we look backwards.
    scanningSpeed *= (endPoint>initialValue ? 1.0 : -1.0);

    // We initialize the return code and other
    // auxiliary variables.
    int returnCode(1), iterations(0);
    double a(initialValue),b(initialValue+scanningSpeed),dx(scanningSpeed);

    // We hunt for a root.
    while (fabs(a)<fabs(endPoint) && iterations < maxIterations)
    {
        // Check for a sign change.
        bool signChange = (function(a)*function(b) < 0.0);

        // If there was a sign change, we use bisection.
        if (signChange)
        {
            double start(a),midPoint,end(b);
            bool asymptote(0);

            while (iterations<maxIterations && !asymptote)
            {
                // Compute the midpiont and the function at that point.
                midPoint = (start+end)/2.0;
                double fMid = function(midPoint);

                // We may have a root.
                if (fMid == 0.0 || fabs(end-start) < toleranceBis)
                {
                    // Check if we have converged on a root or on an asymptote.
                    if (fabs(fMid) < maximumValue)
                    {
                        // We have converged on a root. Use Newton-Raphson.
                        while (iterations<maxIterations)
                        {
                            dx = function(midPoint)/derivative(midPoint);
                            midPoint -= dx;

                            if (fabs(dx) < tolerance)
                            {
                                root = midPoint;
                                returnCode = 0;
                                return returnCode;
                            }
                            else
                                iterations++;
                        }
                    }
                    else
                    {
                        // We have converged on an asymptote.
                        // We jump over it and resume hunting phase.
                        a = b;
                        b = a+scanningSpeed;
                        asymptote=1;
                    }
                }

                // We check the sign at the midpoint. If the function
                // has the same sign at the midpoint and at a end point,
                // we set that end point to the midpoint.
                if (fMid*function(end) < 0.0)
                    start=midPoint;
                else
                    end=midPoint;

                iterations++;
            }
        }
        else
        {
            // We look at the next subinterval.
            a += scanningSpeed;
            b += scanningSpeed;
            iterations++;
        }
    }

    if (returnCode)
    {
        root = sqrt(-1.0);
    }

    return returnCode;
}
